# -------------------------------------------------------------------
# Purpose: Creates Figure B3
# Author:  Max Posch, 25/07/2025
# Usage:   Source this script to generate all panels of the figure
# -------------------------------------------------------------------
# Check that required paths exist
stopifnot(dir.exists(pdataraw))
stopifnot(dir.exists(pdataanalysis))
stopifnot(dir.exists(poutputappendix))

p_unload(tidytable)


## Load data
load(file.path(pdataanalysis, "countyLevel18501940.RData"))
load(file.path(pdataraw, "US_tl2000_cnty_shp_1790_2000.Rdata"))
load(file.path(pdataraw, "US_tl2000_state_shp_1790_2000.Rdata"))

## to address warning about old-style crs
st_crs(shp2000[["1900"]]) <- st_crs(shp2000[["1900"]])
st_crs(state_shp2000[["1900"]]) <- st_crs(state_shp2000[["1900"]])


## simplify shapes
cnty1900 <- shp2000[["1900"]] %>%
        filter(ICPSRSTI != 82, ICPSRSTI != 81, ICPSRSTI != 0) %>%
        select(gisjoin_1900 = GISJOIN, geometry) %>%
        rmapshaper::ms_simplify(keep = 0.05, keep_shapes = TRUE)

state1900 <- state_shp2000[["1900"]] %>%
        filter(ICPSRST != 82, ICPSRST != 81) %>%
        select(geometry) %>%
        rmapshaper::ms_simplify(keep = 0.05, keep_shapes = TRUE)


## Create plots for multiple years
years_to_plot <- c(1850, 1860, 1870, 1880, 1895, 1900, 1905, 1910, 1915, 1920, 1925, 1930)
plot_files <- sapply(years_to_plot, function(yr) {
        cat(paste0("Creating plot for year ", yr, "...\n"))
        plotname <- create_figureB03(yr)
        if (!is.null(plotname)) {
                cat("Figure saved to:", plotname, "\n")
        }
        return(plotname)
})
p_load(tidytable)
